#include <cmath>
#include <limits>
#include <algorithm>
#include <vector>
#include <boost/numeric/interval.hpp>

#include "bounding_box.h"
#include "adaptive_precision.h"
#include "geometry_primitives.h"

using namespace boost::numeric;

namespace 
{
   int sign(const double a) 
   {
      if (a > 0) return  1;
      if (a == 0) return  0;
      else return -1;
   }

   double vect_mul(const segment& s, const point& c)
   {
      return (s.b.x - s.a.x) * (c.y - s.a.y) - (c.x - s.a.x) * (s.b.y - s.a.y); 
   }
   
   double vect_mul_up(const segment& s, const point& c)
   {
      return std::abs((s.b.x - s.a.x) *(c.y - s.a.y)) + std::abs((c.x - s.a.x) *(s.b.y - s.a.y)); 
   }

   bool near_zero(const double mul, const double mul_up)
   {
      double eps = 8 * std::numeric_limits<double>::epsilon() * mul_up;
      return std::abs(mul) <= eps;
   }

   void add_to_adaptive_vector(std::vector<double>& v, double d)
   {
      v.push_back(d);
      std::size_t i = v.size() - 1;
      double res, roundoff;
      while (i > 0)
      {
         res = sum(v[i - 1], v[i], roundoff);
         v[i - 1] = res;
         v[i] = roundoff;
         i--;
      }
   }

   std::vector<double> splitted_adaptive_mul(const segment& s, const point& c)
   {
      double res, roundoff;
      res = mul(s.b.x, c.y, roundoff);
      std::vector<double> v;
      v.push_back(res);
      v.push_back(roundoff);

      res = mul(s.b.x, s.a.y, roundoff);
      add_to_adaptive_vector(v, -res);
      add_to_adaptive_vector(v, -roundoff);
      
      res = mul(s.a.x, c.y, roundoff);
      add_to_adaptive_vector(v, -res);
      add_to_adaptive_vector(v, -roundoff);
      
      res = mul(c.x, s.b.y, roundoff);
      add_to_adaptive_vector(v, -res);
      add_to_adaptive_vector(v, -roundoff);
   
      res = mul(c.x, s.a.y, roundoff);
      add_to_adaptive_vector(v, res);
      add_to_adaptive_vector(v, roundoff);
   
      res = mul(s.a.x, s.b.y, roundoff);
      add_to_adaptive_vector(v, res);
      add_to_adaptive_vector(v, roundoff);
      
      return v;
   }

   int vector_sign(const std::vector<double>& v)
   {
      std:size_t i = 0;
      while (i < v.size()) 
      {
         if (v[i] > 0) 
         {
            return 1;
         } else if (v[i] < 0) 
         {
            return -1;
         }
         i++;
      }
      return 0;
   }

   bool interval_check(const segment& s, const point& c) 
   {
      interval<double> answer_interval = 
         (interval<double>(s.b.x) - interval<double>(s.a.x)) * (interval<double>(c.y) - interval<double>(s.a.y))
         - 
         (interval<double>(c.x) - interval<double>(s.a.x)) * (interval<double>(s.b.y) - interval<double>(s.a.y));
      return zero_in(answer_interval);
   }

   int turn(const segment& s, const point& c)
   {
      double mul = vect_mul(s, c);
      if (near_zero(mul, vect_mul_up(s, c))) 
      {
         
         if (interval_check(s, c)) 
         {
            return vector_sign(splitted_adaptive_mul(s, c));   
         } else 
         {
            return sign(mul);
         }
      } 
      else
      {
         return sign(mul);
      }
   }
}



bool is_bounding_boxes_intersect(const segment& a, const segment& b)
{
   return !((std::max(a.a.x, a.b.x) < std::min(b.a.x, b.b.x) || std::max(b.a.x, b.b.x) < std::min(a.a.x, a.b.x)) ||
                (std::max(a.a.y, a.b.y) < std::min(b.a.y, b.b.y) || std::max(b.a.y, b.b.y) < std::min(a.a.y, a.b.y)));
}



bool are_intersect(const segment& a, const segment& b)
{
   if (is_bounding_boxes_intersect(a, b))
   {
      if (turn(a, b.a) * turn(a, b.b) <= 0 && turn(b, a.a) * turn(b, a.b) <= 0) 
      {
         return true;
      } else
      {
         return false;
      }
   } else
   {
      return false;
   }
}
